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ABSTRACT 

Heat input roughly balances radiative cooling in the gaseous cores of galaxy clusters even when the central 
cooling time is short, implying that cooling triggers a feedback loop that maintains thermal balance. Further¬ 
more, cores with short cooling times tend to have multiphase structure, suggesting that the intracluster medium 
(ICM) becomes locally thermally unstable for cooling times < 1 Gyr. Both observations and theoretical models 
have linked the condensation of cold gas with heating by an active galactic nucleus (AGN) through a cycle in 
which cooling gas fuels the AGN and drives energetic outbursts that reheat the ICM and maintain a state of 
approximate thermal balance. In this work, we use 2D and 3D hydrodynamic simulations to study the onset of 
condensation in idealized galaxy-cluster cores. In particular, we look at how the condensation process depends 
on the ratio of cooling time to freefall time and on the geometry of the gravitational potential. We conclude 
that the ICM can always evolve to a state in which condensation occurs if given enough time, but that an 
initial timescale ratio fcooi/% ^ 10 is needed for thermal instability to grow quickly enough to affect realistic 
cluster cores within a timescale that is relevant for cosmological structure formation. We find that instability 
leads to convection and that perturbations continue to grow while the gas convects. Condensation occurs when 
the timescale ratio in the low-entropy tail of the perturbation distribution drops below fcooi/fif ^ 3, even if the 
volume-averaged timescale ratio is substantially greater. In our simulations, the geometry of the gravitational 
potential does not have a strong effect on thermal stability. Finally, we find that if condensation is powering 
feedback, a conversion efficiency of around 10 “^ for converting the condensed mass into thermal energy is 
sufficient to maintain thermal balance in the ICM. 

1. INTRODUCTION 

X-ray observations of galaxy clusters have revealed that 
the radiative cooling time of gas in many cluster cores 
is much shorter than the Hubble time. If radiative cool¬ 
ing were uncompensated by heating, the gas would radi¬ 
ate away its thermal energy, causing cooling gas to fiow 
toward the center of the cluster. This would be a clas¬ 
sical cooling fiow, in which the accumulating cold gas 
would be observable and would lead to star formation 
rates of > lOOM 0 yr“^ (see Fabian (1994) for a review). 

Instead, X-ray observations reveal little gas cooling be¬ 
low X-ray emitting temperatures (e.g. Peterson et al. 2003; 

Peterson & Fabian 2006) and observed star-formation rates 
that are one or two orders of magnitude lower than pre¬ 
dicted by the classic cooling-fiow model (O’Deaet al. 2010; 

McDonald et al. 2011). Thus, an additional process or 
processes must be heating the ICM to maintain approxi¬ 
mate thermal equilibrium. Several mechanisms have been 
proposed and tested through simulations, including energy 
injection from supernovae (Nagai et al. 2007; Burns etal. 

2008; Skoryetal. 2013), conduction of heat from outside 
of the core (Voigt etal. 2002; Zakamska & Narayan 2003; 

Smith et al. 2013), heating through mergers (Valdarnini 2006; 

Markevitch & Vikhlinin 2007; ZuHone et al. 2010), dynam- 
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ical friction from galaxy cluster motion (Ruszkowski & Oh 
2011; Kim etal. 2005), and feedback from AGN outbursts 
(reviewed by McNamara & Nulsen 2007), which is the mech¬ 
anism we explore in this work. 

AGN feedback is attractive because a simple order-of- 
magnitude estimate shows that an accreting supermassive 
black hole (SMBH) can easily provide enough energy to off¬ 
set cooling. For example, a 10^ Mq SMBH accreting over 
the lifetime of the universe and radiating with a mass-energy 
conversion efficiency of around 10 % would release a total 
of 1062 

ergs, corresponding to an average power output 
of around 10 "^"^ ergs per second—easily enough to offset ra¬ 
diative cooling if a large fraction of that power is injected 
into the ICM (see Churazov et al. (2002) for further discus¬ 
sion). Theoretical and observational studies support the con¬ 
clusion that many cool-core clusters host AGN with enough 
power to balance cooling (e.g., McNamara & Nulsen 2007; 
Dunn & Fabian 2006; Birzan et al. 2004) if a significant frac¬ 
tion of the AGN energy is transfered to the ICM. Nevertheless, 
the details of the AGN fueling process and feedback mode are 
not fully understood. 

If SMBH accretion is to explain thermal regulation 
of the core, then the accretion rate must be linked to 
the thermal properties of the ICM. As pointed out by 
McNamara & Nulsen (2007), if the time-averaged heating 
rate exceeds the cooling rate, the core will heat beyond what 
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is observed, and if it is lower it will fail to prevent gas from 
cooling. More importantly, the short cooling times observed 
in many cluster cores require the heating mechanism to re¬ 
spond on short timescales, on the order of tens of millions 
of years. It is therefore desirable that heating be coupled to 
the cooling rate, to ensure that feedback is able to balance 
cooling both on short timescales and over the lifetime of the 
cluster. Two qualitatively different accretion modes have been 
described in the literature and implemented in numerical sim¬ 
ulations of AGN feedback. Most implementations base the 
black hole accretion rate on the properties of the ambient hot 
gas using modifications of the classic Bondi (1952) analysis 
of smooth, adiabatic accretion, while others rely on conden¬ 
sation and infall of cold clouds to fuel the black hole (e.g., 
Pizzolato & Soker 2005; Gaspari et al. 2012b,a). The analy¬ 
sis of Voit et al. (2014) strongly suggests that the latter “cold 
feedback” mode is more important, because of a universal 
floor observed in the radial cooling-time profiles of galaxy 
clusters that corresponds to the predicted threshold for con¬ 
densation of cold clouds (Sharma et al. 2012a). 

“Cold mode” accretion could be fueled by cold gas con¬ 
densing out of the ICM in response to thermal instability. The 
transition of the ICM from a homogeneous to a hetrogeneous, 
multiphase structure has a long history of investigation us¬ 
ing theoretical arguments and simulations. From a theoretical 
standpoint. Field (1965) studied the evolution of small per¬ 
turbations in cooling plasmas and described an isobaric con¬ 
densation mode, in which variations in temperature and den¬ 
sity may be amplified. Defouw (1970) extended this analysis, 
finding that thermal and convective stability are tightly cou¬ 
pled, a connection further explored in Balbus & Soker (1989). 
The problem of thermal instability in the context of cooling 
flows in clusters was subsequently considered by numerous 
authors (e.g. Cowie et al. 1980; Nulsen 1986; Malagoli et al. 
1987; Loewenstein 1990) who concluded that the cooling 
ICM should indeed be subject to thermal instability. However, 
further analysis by Balbus (1988) and Balbus & Soker (1989) 
using a Lagrangian framework (in contrast to the Eulerian ap¬ 
proach of the earlier works) indicated that the ICM might be 
less susceptible to thermal instability than previously thought, 
especially without the inclusion of a heating term. These 
studies generally take as their starting point an equilibrium 
or steady-state configuration of gas that may not accurately 
capture the behavior of the dynamic ICM. Further, theoretical 
studies are often incapable of dealing with spatially dependent 
heating terms, such as would be expected from star formation 
and AGN feedback. 

There is growing evidence suggesting that the dominant pa¬ 
rameter controlling the transition to a multiphase state and the 
amount of cold gas that condenses is the ratio of gas cooling 
time, tcooh to freefall time, %. Both numerical simulations 
of thermal instability (McCourt et al. 2012) and observations 
of galaxy-cluster cores (Cavagnolo et al. 2008; Rafferty et al. 
2008; Voit et al. 2014) support this conclusion. Without grav¬ 
ity to restore equilibrium, multiphase structure can develop 
within a few cooling times, because collisional cooling pro¬ 
cesses scale with the square of the gas density, allowing 
denser regions to cool faster than their surroundings. If the 
medium is in overall thermal balance, gas clumps that are 
denser than average cool and condense, while underdense re¬ 
gions heat and expand faster than they can cool. In a gravi¬ 
tational potential, however, buoyancy complicates the devel¬ 
opment of thermal instability and can inhibit condensation 
(Malagoli et al. 1987; Balbus & Soker 1989). If the freefall 


time is significantly longer than the radiative cooling time, an 
overdense clump can sink to a denser layer before it can sig¬ 
nificantly cool. 

While theoretical studies provide insight into the general 
physics behind condensation in the ICM, they are necessar¬ 
ily limited by model assumptions and can say little about 
the fate of instabilities that enter the nonlinear regime. In 
recent decades, hydrodynamic simulations such as those of 
Malagoli et al. (1990), McCourt et al. (2012), and Li & Bryan 
(2014) have explored the development of thermal instabil¬ 
ity in astrophysical environments. These works demonstrated 
that condensation can indeed be expected to occur in environ¬ 
ments comparable to the ICM, at a level exceeding the predic¬ 
tions of Balbus & Soker (1989). 

Condensation has been explored in the idealized simula¬ 
tions of McCourt et al. (2012), which show that the growth 
of thermal instabilities is significantly inhibited if 4ooi/% ^ 1- 
However, further studies by Sharma et al. (2012b) have found 
that in a spherical geometry, multiphase gas can still condense 
whenever ^cooi/fif ^10 due to geometric compression (see 
Singh & Sharma (2015) for further discussion.) Gaspari et al. 
(2012b) also finds that a ratio of around 10 is required for 
the formation of cold clumps. Alternately, recent work 
by Li & Bryan (2014) finds that condensation occurs when 
Looi/fif is between 3 and 10. There, condensation is stimu¬ 
lated by interactions between the ICM and an AGN jet. The 
jet entrains cold gas from near the SMBH, pushing it to less 
dense regions. The clump’s positive radial velocity prevents it 
from returning to an equilibrium position, and the gas rapidly 
cools. Finally, observations by Voit & Donahue (2015) find 
that the minimum value of 4ooi/fif in clusters with multiphase 
gas in the form of Ha nebulae generally lies between 5 and 
30. 

In this paper, we use idealized 2D and 3D hydrodynamic 
simulations to study how the onset of condensation depends 
on the ratio of cooling time to freefall time and why there ap¬ 
pears to be a change in cluster core properties around a ratio 
of 10. Section 2 presents simulations based on McCourt et al. 
(2012) in which we explore a wider range of initial condi¬ 
tions. Section 3 analyzes how thermal instabilities grow in 
these simulations and investigate how that growth depends on 
the initial conditions. Section 4 relates this work to previ¬ 
ous theoretical work and discusses the validity of these results 
in the context of real galaxy clusters. Section 5 concludes 
by discussing how future work along these lines may clarify 
how AGN feedback solves the cooling flow problem in galaxy 
clusters. 

2. METHOD 

In this study, we consider simulations of idealized cluster 
cores with planar, cylindrical, and spherical geometries in 2 
and 3 dimensions. The simulations were carried out using the 
AMR Hydrodynamics code Enzo^ (Bryan et al. 2014). Un¬ 
less otherwise noted, 2D runs were conducted on a 300x300 
cell grid with no adaptive mesh, and 3D runs employed a 128^ 
cell root grid with 2 layers of adaptive mesh, with refinement 
based on overdensity, density gradient, and cooling time. We 
do not include magnetic fields, conduction, or the self grav¬ 
ity of the gas. The simulations were analyzed using the yt^ 
analysis toolkit (Turk et al. 201 1). 

^ http://enzo-project.org/ 

^ http://yt-project.org/ 
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2.1. Problem Setup 


We set up the gas in our simulations subject to the constraint 
of hydrostatic equilibrium (HSE) and an Tso-cooling’ initial 
condition, under which the 4ooi/% ratio is uniform throughout 
the volume. Additionally, we run a number of simulations 
using an isothermal initial condition instead of the iso-cooling 
one. 

The setup described in this section applies to all geometries, 
as long as the definition of the height coordinate z changes 
accordingly. In planar geometries, z is the distance from the 
midplane, in cylindrical geometries it is the distance from the 
axis of symmetry, and in spherical geometries it is the distance 
from the origin. We choose a scale height of zs = 100 kpc 
(roughly corresponding to a large cluster), a box size of = 


2z 5', a scale temperature of 7s = 10^ K, and a gravitational 
acceleration scale 


l^l^pZs 


( 1 ) 


so that the gravitational potential energy and thermal energy 
are of similar magnitude at the scale height zs • The cooling 
time is given by 


^COOl(^5 T) 


E 3 nk^T 
|£| 2n^n^k{T) 


( 2 ) 


where E is the thermal energy per unit volume and the form 
of the cooling function A(r) is taken from Sarazin & White 
(1987) for gas of half-solar metallicity. 

The standard normalization of K{T) is used for gas with an 
iso-cooling initial condition of ^cooi/fif = k and we obtain ini¬ 
tial conditions corresponding to other initial values of ^cooi/fif 
by adjusting the normalization of A while keeping the gas 
density and temperature profiles fixed. Two time scales char¬ 
acterizing the initial conditions will be useful in our analysis 
of the onset of thermal instability. One is the cooling time 
Aooi,s at one scale height (z = Zs) at the beginning of the sim¬ 
ulation. The other is the freefall time at one scale height, 
which stays constant throughout the simulation. 

In general, the freefall time of the gas at height z is 


kf{z) = 



( 3 ) 


where the gravitational acceleration defining the potential 
well is 



Z/Zs 


Figure 1. The initial temperature (blue) and density (red) profiles used in this 
work is shown for planar geometry. In simulations with cylindrical and spher¬ 
ical geometries, the gas is isothermal beyond z = 2. The iso-cooling setup, 
with a constant value of fcooi/%5 is shown with a solid line. The isothermal 
setup is shown with a dashed line. Both setups are in hydrostatic equilibrium 
and have a ratio of ^cooi/% = LO at 1 scale height. These profiles are used for 
all simulations; runs with different initial values of ^cooi/% are achieved by 
scaling A(r) after initialization. 


rate at each height. To do this, we sum the total amount of 
cooling in each bin of z, divide by the total volume of the bin, 
and change the sign to get the volumetric heating rate at height 
z. While clearly idealized, this heating prescription ensures 
that the gas remains in overall thermal balance, in agreement 
with the observed thermal behavior of clusters. The validity 
of this prescription is discussed in Section 4.3. 

Eor iso-cooling initial conditions, the initial temperature at 
Zs is 7s • Equations 2 and 3 relate density to temperature via 


hooi ^ 3 nk^T I g(z) 
tff 2A(r)^e^pV 2z 

and the HSE condition for an ideal gas is 


fa 

/imp 




= -piz)g{z) 


( 5 ) 

( 6 ) 


Combining these two expressions gives the temperature 
derivative for iso-cooling initial conditions in the form of an 
ODE: 


giz) = -gs tanh(az/zs) (4) 

for -Rs ^z<Rs and is directed toward either the midplane, 
the symmetry axis, or the origin, depending on the geometry 
of the potential. The relative constancy of g(z) away from the 
origin is meant to mimic the inner region of a spherical gravi¬ 
tational potential in which the mass density is proportional to 
I/z. Eor |z| <^Rs, the tanh function ensures continuity of the 
potential, while the parameter a allows adjustment of its cus- 
piness. In cylindrical and spherical geometries, the magnitude 
of g(z) for z>Rs decreases smoothly decreases smoothly with 
increasing radius to avoid strong acceleration near the bound¬ 
aries. We restrict our analysis to the region -Rs < z < Rs- 
The simulations we present here use a= 1.0, which results in 
a relatively smooth potential with a gradual softening near the 
midplane. 

Eollowing the work of McCourt et al. (2012), we imple¬ 
ment a heating rate that exactly balances the average cooling 


dlnr ^ l nmpg(z)z 1 / ding \1 / dlnA \ 

dlnz [ 2 V dlnz VJvdlnr J 

We integrate this equation to find T (z) and determine the den¬ 
sity from the iso-cooling condition. Eor cylindrical and planar 
simulations, gas outside of Rs is taken to be isothermal and in 
HSE, with r(z) = T{Rs). 

The resulting density and temperature profiles are shown 
in Eigure 1. We impose a temperature fioor at Tfloor = 5.0 x 
10^ K, as we assume that gas below that temperature in¬ 
evitably cools rapidly to much lower temperature. The details 
of the gas fiow below that temperature occur at finer resolu¬ 
tions than are employed in our models, and do not affect the 
overall condensation rate. Einally, we add randomly gener¬ 
ated isobaric perturbations to the gas with an RMS overden¬ 
sity of 0.01 and a fiat spectrum with wave numbers between 2 
and 20 (with k=l corresponding to the box size). The same re¬ 
alization of perturbations is used across all simulations to en- 
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sure consistency. As the gas quickly settles into a convective 
state, the details of the initial perturbations are soon forgotten. 

Figure 1 also shows a comparison between the iso-cooling 
and isothermal setups. In both setups, the initial density and 
temperature profiles do not vary by more than a factor of 4 
throughout the volume. The density is more sharply peaked 
in the isothermal case, leading to a shorter cooling time in 
the center. Consequently, the growth of thermal instabilities 
in the isothermal case is more dependent on height and the 
initial conditions than in the iso-cooling setup. 

3. RESULTS: THE GROWTH OF THERMAL INSTABILITIES 

3.1. Validation of Method 

We begin by conducting simulations with initial conditions 
similar to those in McCourt et al. (2012) to check if our model 
produces qualitatively similar results. Our setup differs from 
theirs in a number of minor details, including the shape of the 
gravitational profile (ours is less cuspy near the midplane) and 
the form of the cooling function A(r). More importantly, the 
region near the midplane does not receive special treatment in 
our simulations, whereas McCourt et al. (2012) shut off heat¬ 
ing and cooling in this region and exclude the midplane re¬ 
gion from analysis. In spite of these differences, we obtain 
qualitatively similar results in the regime 4ooi/fif ^ l-O- We 
see close agreement for the isothermal setup for higher ratios, 
and somewhat more condensation is seen for the iso-cooling 
case for ^cooi/fif ^ l-O- 

Figure 2 shows slices of density in 2D Cartesian simu¬ 
lations with similar initial conditions to those explored by 
McCourt et al. (2012). The top row of Figure 2 uses isother¬ 
mal initial conditions determined by the initial value of ^cooi/fif 
at one scale height. In comparison, the simulations in the 
bottom row use the iso-cooling initial conditions described 
in Section 2 for which ^cooi/% is initially constant throughout 
the simulation volume. The overall behavior is qualitatively 
similar in both cases and resembles the results obtained by 
McCourt et al. (2012) outside of the midplane region. When 
the ratio of timescales is below unity, the gas cools in place 
and forms droplets of condensate that rain down towards the 
midplane. In these cases, convection does not hinder thermal 
instability because the gas is able to adjust its thermal state 
faster than it is able to convect. As the ratio of timescales 
is increased, the dynamics of the gas become increasingly 
dominated by convection, although gas continues to condense 
around the midplane. 

While each vertical pair of models illustrated in Figure 2 
behaves similarly, a number of minor differences can be ob¬ 
served. Principally, condensation occurs more uniformly for 
iso-cooling initial conditions than for isothermal ones. This 
result arises from the differing density profiles needed to sat¬ 
isfy the HSE constraint—isothermal initial conditions have a 
steeper gas-density gradient and consequently a larger range 
in cooling time across the simulation domain. Shorter cool¬ 
ing times near the midplane lead to a ‘cross talk’ effect that is 
more pronounced for isothermal initial conditions. Condensa¬ 
tion of gas near the midplane causes hot, low-density bubbles 
to form there and to rise to greater altitudes, creating inho¬ 
mogeneity at those altitudes on a freefall time scale instead 
of a cooling time scale. This cross talk between lower and 
upper layers complicates the task of interpreting how thermal 
instability and condensation depend on the choice of initial 
timescale ratio at one scale height. The ‘iso-cooling’ con¬ 
dition, while not necessarily more physically valid, reduces 


this cross talk and allows for clearer interpretation of the re¬ 
lationship between the initial timescale ratio and the onset of 
condensation. In contrast, dense midplane gas in models with 
isothermal initial conditions is able to condense quickly even 
when the initial ratio of timescales at one scale height is large. 
This happens because the gas near the midplane will have a 
lower ratio of ^cooi/fif, leading to localized condensation. 

3.2. Instability Growth in the Strong Cooling Regime 

Using the iso-cooling simulations presented in the previous 
section, we have examined the evolution of perturbations for 
the case in which rapid cooling dominates the dynamics of 
the gas. If perturbations are able to cool and collapse more 
rapidly than they can sink, condensation proceeds on a time 
scale ^ ^cooi- When global thermal balance is maintained, 
the average 4ooi/fif of the ambient gas quickly increases as 
condensation lowers the gas mass and density of the ambient 
medium. Within a few cooling times, ^cooi/fif rises to > 10, 
as shown in Figure 3. In this strong-cooling regime, the onset 
of condensation is determined by the growth of the initial per¬ 
turbations and does not depend strongly on the initial ratio of 
Aooi/fif- Condensation continues unabated until the timescale 
ratio is above 10, at which point the cooling is weak enough 
that the condensation rate slows. 

3.3. Instability Growth in the Convective Regime 

When the initial ratio of ^cooi/fif is large, incipient condens¬ 
ing regions sink into the gravitational potential faster than 
they can cool, leading to a roiling, convective state. The con¬ 
vection is subsonic, and although the pressure remains nearly 
constant at a given height, convection does not prevent the 
temperature and density perturbations generated by cooling 
from growing. Figure 4 illustrates the growth of perturbations 
in a medium with an initial timescale ratio of ^cooi/fif = 5.0. 
After 2 cooling times (1 cooling time = 585 Myr), the gas is 
convecting. After 4 cooling times, the gas continues to con¬ 
vect, but the density perturbations have increased. After 6 
cooling times, convection can no longer suppress condensa¬ 
tion of gas near the midplane, and it cools catastrophically. 
After 8 cooling times, a significant amount of the dense gas 
has condensed. 

It is thus clear that the condensation does not simply switch 
on when the average ratio of cooling time to the dynamical 
time drops below some special value. To quantify the transi¬ 
tion of the gas from a relatively smooth, convective state to a 
multiphase medium, we plot in Figure 5 the probability dis¬ 
tribution function of the thermal state of the gas as the run 
with initial ^cooi/fif = 5.0 evolves. After several cooling times 
the distribution of gas in the z-(^cooi/fif) plane has widened 
considerably. After 4 cooling times, gas in the tail of the dis¬ 
tribution has reached a ratio of around 3. At this point, further 
perturbation growth is inevitable and condensation begins. 

Increasing the initial timescale ratio to 4ooi/fif = 20.0 slows 
the condensation process and further restricts it to the mid¬ 
plane region, as shown in Figure 6. Condensation follows the 
same general pattern as in Figure 5, except that it is delayed 
for more than 10 times the initial cooling time and is much 
more pronounced near the midplane. The concentration to¬ 
ward the midplane occurs because cooling gas blobs can set¬ 
tle over a larger number of freefall times and preferentially 
accumulate in the midplane before condensing. 

In all of our simulations, which have iso-cooling initial con¬ 
ditions up to 4ooi/fif = 30, condensation eventually occurs as 
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Figure 2. Slices of gas density are shown for 2D planar simulations with initial values of ^cooi/^ff at one scale height of 0.1, 0.5, 1.0, and 5.0 after the simulation 
has evolved for a time f = 20 fcooi,s • In the top row, the gas is initially isothermal. In the bottom row, the initial timescale ratio is identical throughout the entire 
region. Both models produce qualitatively similar results. In the isothermal case, gas near the midplane has a shorter cooling time than gas above a scale height, 
leading to earlier condensation near the midplane and the creation of hot bubbles that rise up through layers that have not yet begun to condense. In both cases, 
condensation occurs near the midplane in simulations with an an initial value of ^cooi/hf = 5 at one sale height. 



Figure 3. The average ratio of cooling time to freefall time in the ambient 
gas is shown as a function of time for simulations with low initial values of 
fcooi/%- The X axis is in units of the initial cooling time at one scale height, 
^cooi,s 5 rather than absolute time. Values are shown for the 2D planar geometry 
case. Solid lines indicate the volume-averaged value of fcooi/hf within a zone 
0.8 zs <z< L2zs- Dashed lines show the minimum value of tcooi/hf within 
the entire box. At low values of the initial timescale ratio, the gas is able 
to cool in place within a few cooling times, driving the rest of the gas to 
kooi/kf > 10- As the gas cools largely in place, instabilities grow purely on 
the cooling time, leading to similar behavior for all runs. 

long as it is given enough time to develop. Figure 7 shows 
how both the average and minimum values of fcooi/% evolve 
during each run. Condensation in the runs with large val- 



X/Zg X/Zg 


Figure 4. Evolution of gas density in a 2D planar simulation with an iso¬ 
cooling initial condition of koox/kf = 5.0. After two cooling times, the gas is 
clearly convecting. At four cooling times no cold gas has condensed, but the 
amplitude of the perturbations has increased. The perturbations have been 
further amplified after six cooling times, and the first condensate has formed. 
After eight cooling times, the densest gas near the center has entered into 
runaway cooling, leading to continuous condensation. 

ues of tcoo\/tfi may be surprising in light of recent theoreti¬ 
cal studies predicting that the medium should become multi¬ 
phase only if tcoo\/tfi <10 (Sharma et al. 2012b; Gaspari et al. 
2012b; Singh & Sharma 2015), and we will discuss possible 
explanations for this difference in Section 4. 

Figure 8 shows the same simulations as Figure 7 plotted 
with time in units of the freefall time at one scale height rather 
than the initial cooling time. As all simulations use the same 
gravitational potential, the freefall time is a standard clock 
and corresponds to the same time interval in physical units 


Density (g cm' 
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Figure 5. Evolution of the mass-weighted probability distribution for the 
ratio of cooling time to freefall time for a 2D planar geometry with initial 
fcooi/hf = 5.0. The dashed black line shows the volume-weighted average ra¬ 
tio as a function of height. Note that when gas condenses, most of the volume 
is occupied by the hot gas, meaning that the volume-averaged ratio will tend 
to lie above the mass-weighted mean. The first panel shows the initial state 
of the gas where the timescale ratio is held constant throughout (with some 
spread due to the initial perturbation spectrum). Att = 2.0 tcoo\,s^ the gas has 
entered into a convective state and although condensation has not yet com¬ 
menced, a spread in gas properties is evident. By t = 6.0 ^cooi,s, a portion of 
the gas has reached a state with fcooi/hf ~ 2-3, and the condensation process 
has begun. Although some gas is entering into the cold phase, the volume av¬ 
eraged ratio of fcooi/% remains near its initial value as the cold gas occupies 
negligible volume. 



Figure 7. Same as Figure 3, except for runs with larger initial values of 
^cooi/%- In sach simulation, the minimum value of fcooi/% decreases on a 
timescale roughly proportional to the cooling time. When the initial ratio 
is higher, it takes several cooling times for gas to develop regions with a 
minimum timescale ratio near unity; therefore, condensation is delayed in 
these runs. Note that runs with larger values of fcooi/hf have a low overall 
cooling rate which, combined with the temperature floor of Tfloor = 5x10^ 
K, produces the floor in the timescale ratio. 



Figure 6. Same as Figure 5 except for an initial timescale ratio of ^cooi/^ff = 
20.0. By r = 5.0 kooh perturbations have started to grow but have not yet led 
to condensation. As the gas is able to undergo more freefall times per cooling 
time than in the case of kooi/hf = 5.0, cooler gas is able to effectively settle 
towards the midplane. Nevertheless, condensation is still able to occur near 
the midplane even though the volume averaged value of kooi/hf remains near 
10 . 

(approximately 117 million years). When the initial ratio of 
cooling time to freefall time exceeds ^10, condensation oc¬ 
curs after ^100 freefall times, corresponding to a timescale 
comparable to the Hubble time. 

3.4. Transition to the Condensed State 

While studying the growth of thermal instabilities gives in¬ 
sight into the conditions under which gas will condense, it 
does not necessarily explain how the gas reaches a condensed 
state or how an individual parcel of gas behaves. Figure 9 de¬ 



Figure 8. Same as Figure 7 except with time plotted in units of the freefall 
time at one scale height rather than the initial cooling time at one scale height. 


picts the gas distribution in the p - T plane integrated over 30 
cooling times. To compute this distribution, we bin gas mass 
m p-T space in each data output (which are evenly spaced 
in time), sum over all of the outputs, and normalize so that 
the integral over the distribution is equal to 1. This probabil¬ 
ity distribution corresponds to the probability of a parcel of 
gas being found in a given thermodynamic state at some point 
during the simulation. The gas is for the most part constrained 
to a line of constant pressure with spread due to gravitational 
stratification. The distribution has two peaks; a low density. 































































Growth and Evolution of Thermal Instabilities 


7 





Figure 9. The probability distribution function of the gas in the p-T plane, 
averaged over 30 cooling times in the 2D planar simulation with an initial 
timescale ratio of fcooi/% = 5. Note that a large fraction of the gas is located 
at the temperature floor, near the x-axis at a density slightly above 10“^^ g 
cm“^. As the convection and condensation processes proceed subsonically, 
the process is largely isobaric, with a modest spread due to gravitational strat- 
iflcation. Lines of constant cooling time are shown as dashed lines, labeled 
with the ratio of cooling time to freefall time at 1 scale height. Note that the 
gas spends very little time between the line tcooi/kf = 2 and the temperature 
floor, indicating that once the threshold is reached, condensation proceeds 
rapidly. 

high temperature node in which the gas is convecting, and a 
cool, low temperature node representing the condensed state 
of the gas. The probability of finding gas in the connecting 
region is low, indicating that condensation from the hot phase 
into the cold phase proceeds rapidly once it begins. 

Figure 10 illustrates the dynamics of the gas during the con¬ 
vective stage and the condensation process using the motion 
of a Lagrangian tracer particle in phase space. The figure 
shows the path of a representative particle which condenses 
early in the simulation. For several cooling times, the gas 
simply convects within a narrow portion of phase space. As 
the thermal perturbations are amplified, the gas is driven to 
a colder, denser state which is where condensation occurs. 
When the gas does condense, the condensation process is very 
rapid, and the gas stays in the condensed phase afterwards. 


3.5. Condensation Rate 

In our simulations, condensed gas remains in the condensed 
state and settles towards the center. After the onset of conden¬ 
sation the gas segregates into two phases - the cool condensed 
material in the center and the hot, convective gas that remains 
uncondensed. This departure from the expectation of self¬ 
regulation is a consequence of our feedback implementation 
and is discussed further in Section 4.4. Still, it is instructive 
to examine the rate of condensation in our simulations, as is 
shown in Figure 11 . Following the onset of condensation, in¬ 
stabilities continue to grow on the cooling timescale. Each 
simulation behaves similarly on a thermal timescale, with a 
roughly linear growth in the total condensed fraction. 




(0 

Q. 

E 


10 ^ 




t/fool 


Figure 10. The dynamics of fluid during the condensation process are shown 
in the dynamics of a Lagrangian tracer particle through phase space. The 
particles are inserted during initialization in the 2D planar simulation with 
an initial timescale ratio of 5. The upper left panel shows the particle’s path 
through p-T space, with the color of the line showing elapsed time in coohng 
times. The upper right panel also shows the path through p-T space, but is 
colored by the ratio of cooling to freefall time. The bottom left panel shows 
particle height vs. time, and the bottom left shows the timescale ratio as a 
function of time. 



Figure 11. The fraction of mass in the condensed state is shown as a function 
of time for 2D planar runs with large initial values of ^cooi/^ff- The condensed 
fraction is measured over the entire domain. 


3.6. Effect of Geometry 

Simulations with different geometries are shown in Figure 
12. All simulations use the same initial conditions (fcooi/fif = 
5.0). In the spherical case, gravity pulls towards the origin, 
while in the cylindrical setups gravity pulls towards the sym¬ 
metry axis. All simulations exhibit similar thermal behavior. 


After 10 cooling times, the gas has entered into a convective 
state and condensation has begun near the center of the poten¬ 
tial. In the non-planar runs, less gas condenses as the region 
near the center occupies less volume. Nevertheless, we do 
not observe a significant change in the condensation process 
among simulations with different geometries. 
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Density [g cm^ ] 

10 10 10 



- 2-10 1 - 2-10 1 2 

X/Zg X/Zg 

Figure 12. The evolution of runs with an initial fcooi/hf of 5.0 are shown 
for different geometries at r = lOfcooi- All runs use an identical setup with 
respect to the radial coordinate, though the dehnition or the radial coordinate 
is changed based on the geometry of the simulation. The 2D runs use a static 
grid of 300 x 300 cells, while the 3D runs use a root grid of 128^ cells with 
2 layers of adaptive mesh. The slices of density through the cylindrical runs 
are taken perpendicular to the axis of symmetry. 

4. DISCUSSION AND RELATIONSHIP TO RELATED WORK 

Our simulations would seem to indicate that any medium 
subject to a heating/cooling balance as we have described in 
our model will eventually succumb to thermal instability and 
produce condensation. Nevertheless, observations seem to in¬ 
dicate that clusters with time scale ratios above roughly 10 do 
not produce much multiphase gas. To explain this discrep¬ 
ancy, we note that at a radius of around 30 kpc, a time scale 
ratio of 10 in a large galaxy cluster corresponds to a cool¬ 
ing time on the order of a Gyr. Physical processes such as 
mergers, star formation, and AGN feedback occur on shorter 
timescales, rendering the condensation process sub-dominant 
in these cases. Therefore, in a realistic cluster environment 
only clusters with a cooling time to freefall time ratio of < 10 
are likely to develop condensation. An important caveat to 
this observation is that while the growth of thermal instabili¬ 
ties from initially small perturbations may be unimportant on 
cluster timescales, if the gas is inhomogeneous due to other 
physical processes (such as an AGN jet) condensation may 
occur in the tail of the thermal distributions shown in Figure 
5 and 6. Thus, predicting the onset of condensation is not as 
simple as measuring the value of fcooi/%; the level of inhomo¬ 
geneity must also be taken into account. 

Our simulations examine the formation of multiphase gas in 
an idealized setting wherein global balance between heating 
and cooling is strictly enforced. While this model gives rise 
to results that are qualitatively consistent with observations, 
it clearly neglects the complex physics of AGN feedback and 
heat transport which occur in real clusters. In this section, we 
discuss our results in light of current observations of multi¬ 
phase gas and previous simulations of condensation and con¬ 
sider the complications that inclusion of additional physical 
processes would cause. 

4.1. Observations of Multiphase Gas 


Owing to the timescales involved and the limits of current 
telescopes, astronomers can not directly observe the conden¬ 
sation process in the ICM. Nevertheless, the past decade has 
deepened the field’s appreciation of a fascinating dichotomy 
in cluster properties when cluster cores are probed for cold 
gas and signatures of AGN feedback. While cooling is gen¬ 
erally suppressed in cool-core galaxy clusters (Peterson et al. 
2003; Peterson & Fabian 2006), at least some cold gas is ob¬ 
served in galaxies with low central temperatures, as seen in 
the works of McDonald et al. (2010) and Werner et al. (2010). 
Cavagnolo et al. (2008) considers the entropy profiles, radio 
emissions, and presence of Ha in the ACCEPT sample of 222 
galaxy clusters. As Ha emission requires the presence of cold 
(relative to the ICM) gas, the presence or absence of Ha in a 
cluster may be taken as an indicator of multiphase gas. In the 
clusters with Ha observations, Ha is conclusively detected 
in slightly over half of the sample. A strong correlation is 
seen between the presence of Ha and the core entropy; clus¬ 
ters with Ha have central entropies below 30 keV cm^, while 
those without Ha detections tend to lie above the 30 KeV line. 
When the entropy profile is used to infer a cooling time, (as 
in Voit & Donahue 2015; Voit et al. 2014) a central entropy 
of 30 KeV corresponds to a cooling time of around 1 billion 
years, consistent with a cooling time to freefall time ratio of 
around 10. In clusters in the ACCEPT sample, those with 
detected Ha emission consistently have 4ooi/% values below 

20, while those without Ha detections lie entirely above 
that value. 

4.2. Simulations of Multiphase Gas 

McCourt et al. (2012), upon which this study is based, finds 
that precipitation will occur rapidly if the gas is able to cool in 
place, which occurs when ^cooi/fif ^ 1- The authors also con¬ 
clude that the condensation process is relatively insensitive to 
variations in the heating rate and mechanism. Employing a 
similar method, Sharma et al. (2012b) finds that condensation 
may occur in gas with a timescale ratio of < 10 in a spherical 
simulation, an enhancement they attribute to the compression 
of overdense blobs descending in a spherical geometry. Ad¬ 
ditionally, Sharma et al. (2012b) concludes that condensation 
does not occur when the timescale ratio rises above 10. 

Analytic work has lent further credence to the idea that 
^cooi/% ^10 represents a critical threshold for condensa¬ 
tion. Singh & Sharma (2015), extending the analysis of 
Pizzolato & Soker (2005), finds that small instabilities may 
grow when 4ooi/fif ^ 1 for planar geometries and, when the 
effects of geometric compression are included, may grow for 
^cooi/fif ^ 10 for spherical geometries. While the results pre¬ 
sented in Section 3 of this paper suggest a moderately higher 
threshold for the planar case, we believe that these results are 
largely consistent with Singh & Sharma (2015) in the context 
of individual overdensities cooling and condensing in place. 
In our simulations, however, we see that overdensities in a 
medium above the critical threshold oscillate, leading to a 
roiling state that develops further perturbations. This cross 
talk effect between layers generates non-linear perturbations 
and causes the temperature dispersion in the medium to grow 
on the cooling timescale. When the cold tail of the distribu¬ 
tion has dropped to ^cooi/fif ^2-3 the condensation process 
begins. Thus we find that the mechanism responsible for con¬ 
densation above the critical threshold is not geometric com¬ 
pression but the continued growth of perturbations following 
the onset of convection in the gas. 

Simulations that employ more realistic heating mechanisms 
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also find that condensation occurs in galaxy clusters, albeit 
under somewhat different circumstances than in simulations 
with idealized heating. Li & Bryan (2014) employ an AGN 
feedback algorithm in which heating is triggered by cold gas 
accretion. The study finds that condensation occurs when 3 < 
^cooi/fif ^10- This condensation occurs along the axis of the 
jet, where dense gas is dragged up and is able to cool as it 
falls. This is consistent with our findings, in which the thermal 
instability can grow when ^cooi/% ^10. but only when gas 
is sufficiently hetrogeneous. Similarly, Gaspari et al. (2012b) 
employs jet heating in response to accretion and finds that 
multiphase gas can form when ^cooi/fif ^10- 

4.3. Caveats and Limitations 

In this study, we have used an idealized model to simulate 
the onset of condensation in galaxy clusters. While the sim¬ 
plicity makes this model easy to analyze, we have left out 
physics that may have significant impact on the development 
of a multiphase medium. In particular, conduction and the 
presence of magnetic fields may inhibit or shape the growth of 
condensation. Conduction works to smooth out temperature 
perturbations, while magnetic fields will lead to conduction 
being anisotropic. 

Magnetic fields in clusters are poorly understood. While 
weak, they are known to be present and may be dynamically 
important in cluster cores (Carilli & Taylor (2002) and refer¬ 
ences therein). More important for this work, magnetic fields 
in a plasma will lead to anisotropic conduction, channeling 
heat along the direction of magnetic field lines as explored in 
Ruszkowski et al. (2011). Similarly, Wagh et al. (2014) stud¬ 
ies the growth of thermal instabilities in a spherical setup and 
includes both conduction and magnetic fields. Anisotropic 
conduction is not found to inhibit condensation, but does lead 
to the formation of filaments rather than globules of dense gas. 
Conduction is found to inhibit condensation if the efficiency 
is above 0.3 of the Spitzer value. 

In addition to these omissions, our assumed heating func¬ 
tion does not capture the true physical process responsi¬ 
ble for transferring energy from the AGN to the ICM. 
While the details of AGN heat transfer are not currently un¬ 
derstood, several mechanisms have been proposed, includ¬ 
ing shocks (McNamara & Nulsen 2012; Ruszkowski et al. 
2004), cosmic rays (Sharmaetal. 2010; Fujitaetal. 2013; 
Fujita & Ohira 2011, 2012), turbulent mixing (Sharma et al. 
2009; Banerjee & Sharma 2014), PdV work from the infia- 
tion of hot bubbles (McNamara & Nulsen 2007; Birzan et al. 
2004), and the uplifting of cool gas by rising bubbles 
(Million et al. 2010). The actual heating function is unlikely 
to maintain perfect thermal balance, and presumably does not 
act in a strictly volumetric sense as assumed in this work. 
Still, the lack of cold gas and star formation in cool-core 
clusters implies that the heating function must broadly main¬ 
tain thermal equilibrium, making the model considered in this 
work physically relevant. 

4.4. Self-Regulation 

The gas does not reach a steady state, as might be expected 
for an ideal self-regulating system. Instead, condensation 
continues in the convective gas after the condensation process 
has begun, increasing the separation between the hot and cold 
phases. In real clusters, feedback is expected to operate in a 
thermostat-like manner, which should produce a rough ther¬ 
mal equilibrium. The lack of self-regulation in our simula¬ 
tions is purely an effect of the heating model that we employ. 



Figure 13. The feedback efficiency necessary to maintain thermal balance 
in the hot gas is shown for 2D planar runs starting with different values of 
tcooi/%- The required efficiency is calculated as e = E/McoidC^^ where E is 
the cooling rate of all gas above the temperature floor. Both the cooling rate 
and the condensation rate have been smoothed over a cooling time. 

and does not accurately capture the response of feedback to 
condensation. However, if we imagine feedback to be pow¬ 
ered by condensation, we can use the calculated heating rate 
to determine what feedback efficiency would be necessary for 
the system to balance radiative losses. 

During accretion, AGN are expected to convert a significant 
fraction of the infalling mass into energy that is then returned 
to the surrounding medium. The feedback rate can be related 
to the mass accretion via 

E « eMc^ (8) 

where E is the total energy output, e is an efficiency parame¬ 
ter, and M is the mass accretion rate. Under the assumptions 
that all of the condensing gas is used to power feedback and 
that all of the feedback energy is transferred to the ICM, we 
have estimated the conversion efficiency necessary to main¬ 
tain thermal balance. The estimate is shown in Figure 13 for 
several initial values of fcooi/fif • Once condensation has begun, 
the required efficiency in all runs reaches a value of around 
10“^, in line with the values found in Sharma et al. (2012b). 
As the accumulation of cold gas near the midplane is an ar¬ 
tifact of our setup and would not be expected in a real clus¬ 
ter, we calculate the cooling rate over the ambient gas, which 
is that gas that is above the temperature fioor of 5.0 x 10^ 
K. Although we do not explore the mechanism for releasing 
mass-energy from condensed gas in this work, if condensa¬ 
tion resulting from the growth of thermal perturbations is in 
principle capable of balancing radiative cooling in the ICM, 
thermal instability must be taken seriously as a feature of a 
self-regulating energy cycle in cool-core clusters. 

5. CONCLUSIONS AND FUTURE WORK 

In this study, we have investigated the onset of convection 
in a thermally unstable medium using an idealized model, 
including a heating scheme that strictly enforces a global 
heating-cooling balance. Although a simplification, this 
model gives insight into the conditions necessary for the onset 
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of condensation in a gravitationally stratified medium such as 
that in a cool-core galaxy cluster. This study indicates that 
condensation proceeds as follows: 

• If heating is able to balance cooling at all radii, thermal 
instabilities will grow in amplitude, regardless of the 
initial conditions. 

• If the ratio of the cooling to the freefall time is < 2, (the 
strong cooling regime) the gas will condense in place, 
driving the volume-averaged ^cooi/fif value above 10. 

• Above a ratio of ^cooi/fif ~ 10, perturbations will grow 
on a timescale proportional to the cooling time. 

• Once the perturbation distribution has broadened, gas 
with ^cooi/fif ~ 2-4 will condense, even if the volume- 
averaged ratio of 4ooi/fif is above 10. 

• If the timescale ratio is > 10, the timescale for conden¬ 
sation to occur in gas with ^cooi ^ 1 Gyr is comparable 
to the Hubble Time and greatly exceeds other relevant 
cluster timescales. 

A fundamental limitation of this work is that the model as¬ 
sumes a heating function that is idealized and does not mimic 
a specific physical process. In preparation for future work, 
it will be necessary to examine a greater variety of heating 
modes, including models more analogous to jet feedback and 
quasar winds from accreting supermassive black holes. The 
physical processes underlying black hole accretion, feedback, 
and heat transfer to the ICM are still poorly understood, and 
elucidating them will form the focus of future studies. 
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